#-------------------------------------------------------#
#   Quality Assessment of the Academic Freedom Index    #
#-------------------------------------------------------#

# Title: "Quality Assessment of the Academic Freedom Index: Strengths, Weaknesses, and How Best to Use It" 
# Authors: "Lott, Lars", "Spannagel, Janika"
# date: 2024-09-19
# journal: Perspectives on Politics
# DOI: TBA
# written under "R and RStudio (4.4.1)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# load packages
library(here)
library(tidyverse)
library(ggpubr)
library(countrycode)
library(ggrepel)
library(ggthemes)
library(lubridate)
library(dotwhisker)
library(lavaan)

library(coda)
library(MASS)
library(rjags)
library(parallel)
library(magrittr)
library(dplyr)
library(tools)
library(devtools)
library(tinytable)

library(lavaan)
library(stargazer)
library(vutils)

set.seed(1234)


#### Traditional CFA (Section 3 in Appendix)

vnames <- c("v2cafres", "v2cafexch", "v2cainsaut", "v2casurv", "v2clacfree")

#### Load V-Dem dataset for frequentist factor analysis ##

vdem_full_v13 <- readRDS(file.path("data", "vdem_v13", "vdem_13_cy", "V-Dem-CY-Full+Others-v13.rds"))


## Fit model 1-dim model
dim1 <- cfa('f =~ v2cafres + v2cafexch + v2cainsaut + v2casurv + v2clacfree',
            data=vdem_full_v13, std.lv=TRUE, missing="ML")

## Loadings are Std.all column for f, Uniquenesses are Std.all columnv for Variances
## Reports TLI, RMSEA, and SRMR
summary(dim1, fit.measures=TRUE, standardized=TRUE)


## First fit 1-dim model without imputation, to compare
dim1lwd <- cfa('f =~ v2cafres + v2cafexch + v2cainsaut + v2casurv + v2clacfree',
               data=vdem_full_v13, std.lv=TRUE)
summary(dim1lwd, fit.measures=TRUE, standardized=TRUE)

### Two dim models

dim2 <- cfa('f1 =~ v2cafres + v2cafexch + v2clacfree 
             f2 =~ v2cainsaut + v2casurv',
            data=vdem_full_v13, std.lv=TRUE)
summary(dim2, fit.measures=TRUE, standardized=TRUE)



## Test null of identical model fit
anova(dim2, dim1lwd) # summarized begin last para appendix section 3



## Extract Loadings and Uniquenesses ##

dim1_container <- lavInspect(dim1lwd, "Std.all")

loadings_dim1 <- dim1_container$lambda
uniquenesses_dim1 <- diag(dim1_container$theta)


dim1_table  <- tibble(data.frame(Measure = c("v2cafres", "v2cafexch", "v2cainsaut", "v2casurv", "v2clacfree"), 
                          Loadings = loadings_dim1, 
                          Uniqueness = uniquenesses_dim1))
dim1_table$f <- round(dim1_table$f, 3)
dim1_table$Uniqueness <- round(dim1_table$Uniqueness, 3)

tinytable::tt(dim1_table) %>% save_tt("outputs_original_data/Table_B1.tex", overwrite = T)

## Extract Loadings and Uniquenesses ##

dim2_container <- lavInspect(dim2, "Std.all")

loadings_dim2 <- dim2_container$lambda
uniquenesses_dim2 <- diag(dim2_container$theta)


dim2_table  <- tibble(data.frame(Factor = c("Dimension 1", "", "", "Dimension 2", ""),
                                 Measure = c("v2cafres", "v2cafexch", "v2clacfree", "v2cainsaut", "v2casurv"), 
                                 Loadings = loadings_dim1, 
                                 Uniqueness = uniquenesses_dim1))
dim2_table$f <- round(dim2_table$f, 3)
dim2_table$Uniqueness <- round(dim2_table$Uniqueness, 3)

tinytable::tt(dim2_table) %>% save_tt("outputs_original_data/Table_B2.tex",overwrite = T)



